OLD SHAPE¶
Apenas eco-regiões
In [999]:
%%HTML
<h1> Índice</h1>
<h3>
<a href="#1">clicar aqui!</a> para ir para as primeiras figuras
</h3><br>
<h3> outras figuras:</h3><br>
<a href="#varintermonbio">Variação inter-anual por mês em cada bioma</a>
<br>
<a href="#varintrabio"> Variação intra_anual em cada bioma</a>
In [1]:
# imports
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as cls
from pandas import Series, DataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset
import matplotlib.ticker as ticker
from collections import defaultdict, namedtuple
from pandas import ExcelWriter
import xlsxwriter
import pandas as pd
In [694]:
import calendar
In [474]:
plt.rcParams["font.family"] = "Times New Roman"
In [3]:
def get_ecos():
m = Basemap(projection='cyl')
shpf = m.readshapefile('ecoregions2017/ecoregions2017','ecos',linewidth=0.1)
return m.ecos_info, m.ecos
In [4]:
months = "Jan Fev Mar Abr Mai Jun Jul Ago Set Out Nov Dez".split(' ')
ecos_info, ecos = get_ecos()
Eco = namedtuple('Eco', ['name', 'biome_name', 'biome_num', 'area'])
info_dict = {eco['ECO_ID']:Eco(eco['ECO_NAME'],eco['BIOME_NAME'],eco['BIOME_NUM'], eco['SHAPE_AREA']) for eco in ecos_info}
eco_names = []
for i in range(847):
name = info_dict[i].name
if type(name) == bytes:
name = name.decode('latin-1')
eco_names.append(name)
biome_names = []
for i in range(847):
biome_names.append(info_dict[i].biome_name)
biome_nums = []
for i in range(847):
biome_nums.append(info_dict[i].biome_num)
areas = []
for i in range(847):
areas.append(info_dict[i].area)
In [9]:
biome_nums[0] = 0.0
In [5]:
result = {}
for i in range(1,18):
suf = '20{0:02d}'.format(i)
in_file= 'result' + suf
result[suf] = DataFrame(np.load(in_file), columns =months)
In [6]:
for year in result:
result[year].index.name ='Eco_Id'
In [11]:
for year in result:
result[year] = result[year].transpose()
In [14]:
head = DataFrame()
head = head.assign(ECO_NAME = eco_names, BIOME_NUM=biome_nums,BIOME_NAME=biome_names,AREAS= areas).transpose()
In [15]:
result['Atributos'] = head
In [17]:
df_finall = pd.concat(result)
In [18]:
df_finall
Out[18]:
Não dividido por área¶
Total por Bioma¶
In [25]:
df_finall.loc['Atributos'].loc['BIOME_NUM']
Out[25]:
In [32]:
df_finall[:][:12*17]
In [33]:
biome_grouped = df_finall[:][:12*17].groupby(df_finall.loc['Atributos'].loc['BIOME_NUM'], axis=1)
In [34]:
biome_grouped.sum()
Out[34]:
In [35]:
#writer = ExcelWriter('results_newstyle2.xlsx',engine='xlsxwriter')
#df_finall.to_excel(writer,sheet_name='Burned Pixels',engine='xlsxwriter')
#writer.save()
In [37]:
total_biome = biome_grouped.sum().sum()
RESULTADO¶
In [38]:
total_biome
Out[38]:
Total por Ecoregião¶
In [41]:
total_eco = df_finall.iloc[:12*17].sum()
In [42]:
total_eco
Out[42]:
Dividido por área¶
Total por ecoregião¶
In [49]:
total_eco_area = total_eco/(areas)
Total por Bioma¶
In [53]:
df_finall[:][12*17:-1:2]
Out[53]:
In [57]:
df_finall[:][12*17 +2: 12*17 +3]
Out[57]:
In [63]:
df_finall.loc['Atributos'].loc['AREAS']
Out[63]:
In [64]:
df_finall.loc['Atributos'].loc['AREAS'].groupby(biome_nums).sum()
Out[64]:
In [82]:
df_finall.T[df_finall.T['Atributos']['BIOME_NUM'] == 14]['Atributos']['AREAS'].sum()
Out[82]:
In [83]:
areas_biome = df_finall.loc['Atributos'].loc['AREAS'].groupby(biome_nums).sum()
In [89]:
areas_biome
Out[89]:
RESULTADO¶
In [91]:
total_biome_area = total_biome/areas_biome
In [88]:
total_biome.iloc[14]/areas_biome[14]
Out[88]:
In [90]:
302628.0/27.665056
Out[90]:
SOMAS MENSAIS¶
Por ecoregião¶
In [119]:
df_finall[1: 12*17:12]
Out[119]:
In [113]:
months
Out[113]:
In [137]:
monthly = [ df_finall[n: 12*17:12].sum()/areas for n in range(12)]
RESULTADO¶
In [138]:
monthly
Out[138]:
In [139]:
dict_monthly = dict(zip(months,monthly))
In [140]:
dict_monthly['Jan']
Out[140]:
Por Bioma¶
In [149]:
monthly_biome = [ df_finall[n: 12*17:12].sum() for n in range(12)]
In [150]:
monthly_biome = [mb.groupby(biome_nums).sum()/areas_biome for mb in monthly_biome]
RESULTADO¶
In [152]:
monthly_biome
Out[152]:
In [153]:
dict_monthly_biome = dict(zip(months,monthly_biome))
SOMAS ANUAIS¶
Por Ecoregião¶
In [154]:
annual = [ df_finall[n*12: (n+1) *12].sum()/areas for n in range(17)]
In [155]:
annual
Out[155]:
In [158]:
years = [str(y) for y in range(2001,2018)]
In [159]:
years
Out[159]:
In [160]:
dict_annual = dict(zip(years,annual))
In [161]:
dict_annual['2002']
Out[161]:
Por bioma¶
In [162]:
annual_biome = [ df_finall[n*12: (n+1) *12].sum() for n in range(17)]
In [165]:
annual_biome = [ ab.groupby(biome_nums).sum()/areas_biome for ab in annual_biome]
In [167]:
dict_annual_biome = dict(zip(years,annual_biome))
Gráficos totais por ano e por meses¶
sem áreas¶
In [829]:
annual_total = [ df_finall[n*12: (n+1) *12].sum().sum() for n in range(17)]
fig, ax = plt.subplots()
ax.plot(years, annual_total)
labels = ax.get_xticklabels()
plt.tight_layout()
plt.setp(labels, rotation=30, fontsize=10);plt.show()
In [182]:
monthly_total = [ df_finall[n: 12*17:12].sum().sum() for n in range(12)]
In [223]:
fig, ax = plt.subplots()
#ax[0].plot( np.arange(1,13),monthly_total)
#labels = ax.get_xticklabels()
#plt.setp(labels, rotation=30, fontsize=10)
#plt.xticks(months,months)
ax.plot( np.arange(12),monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months)
plt.show()
In [185]:
months
Out[185]:
In [230]:
all_months = [df_finall.iloc[m].sum() for m in range(17*12)]
In [237]:
import itertools
ym = itertools.product(years,months)
ymlabels = [y + m for y, m in ym]
ymlabels
from statistics import mean
media = mean(all_months)
fig, ax = plt.subplots()
ax.plot( np.arange(12*17),all_months)
ax.axhline(y=media, color='r', linestyle='--')
ax.set_xticks(np.arange(12*17))
ax.set_xticklabels(ymlabels)
plt.setp(labels, rotation=60);
plt.show()
Out[237]:
In [264]:
df_finall
Out[264]:
In [265]:
annual_total
Out[265]:
In [281]:
n = 0
In [282]:
df_finall[n*12: (n+1) *12].sum().sum()
Out[282]:
In [288]:
ar_annual_biome = np.array(annual_biome)
In [592]:
#testing
lixo = np.array(annual_biome)
In [593]:
lixo.shape
Out[593]:
In [293]:
ar_annual_biome.reshape(17,15)
Out[293]:
In [304]:
anbio = ar_annual_biome.T
In [305]:
anbio[0]
Out[305]:
In [301]:
annual_biome
Out[301]:
In [318]:
fig, ax = plt.subplots(5,3, figsize =(18,12))
for i in range(5):
for j in range(3):
ax[i][j].plot(years, anbio[i*3 + j])
labels = ax[i][j].get_xticklabels()
plt.setp(labels, rotation=30, fontsize=10)
plt.show()
In [319]:
basket = 'BasketSemNomes.csv'
In [327]:
basket_df = pd.read_csv(basket,sep=';')
In [330]:
basket_df
Out[330]:
In [332]:
basket_df.describe()
Out[332]:
In [331]:
!cat ITDcluster.r
TODOS OS PIXEIS¶
In [334]:
allpixels = np.load('../MODIS/MCD64A1/allpixels.npy')
In [337]:
allpixels[0]
Out[337]:
In [338]:
annual_total[0]
Out[338]:
In [340]:
an_total = np.array(annual_total[:-1])
In [344]:
diffs = (allpixels - an_total)*100/allpixels
In [350]:
diffs.max()
Out[350]:
In [349]:
fig, ax = plt.subplots(3,1)
ax[0].plot(years[:-1], annual_total[:-1])
ax[0].set_title("Pixeis havendo selecção p/ ecoregião")
ax[1].plot(years[:-1], allpixels)
ax[1].set_title("Pixeis não havendo selecção p/ ecoregião")
ax[2].plot(years[:-1], diffs)
ax[2].set_title("Diferença em percentagem %")
labels = [ax[i].get_xticklabels() for i in range(3)]
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
biomes again¶
In [357]:
biome_names_dict = dict( (num, name) for num, name in zip(biome_nums, biome_names))
In [358]:
biome_names_dict
Out[358]:
In [364]:
fig, ax = plt.subplots(5,3, figsize =(18,12))
for i in range(5):
for j in range(3):
ax[i][j].plot(years, anbio[i*3 + j])
labels = ax[i][j].get_xticklabels()
ax[i][j].set_title(biome_names_dict[i*3 + j])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
todos os gráficos com a mesma escala?¶
experimentar aplicar a escala das savanas às florestas mediterrânicas
In [365]:
savanas = 7; mediterraneo = 12
In [368]:
fig, ax = plt.subplots(2,1, sharey=True)
ax[0].plot(years, anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(2,1, sharey=True)
ver também sharey='col' e sharey='row'
In [375]:
log_anbio = (np.log(anbio) - np.log(anbio.min()))/(np.log(anbio.max()) - np.log(anbio.min()))
In [376]:
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escala logarítmica normalizada')
ax[0].plot(years, log_anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, log_anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
In [380]:
log2_anbio = np.log(anbio -anbio.min() + 1) / np.log(anbio.max() - anbio.min() +1 )
In [381]:
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escala logarítmica normalizada 2')
ax[0].plot(years, log2_anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, log2_anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
Comparar fórmulas de logaritmos¶
In [382]:
x = np.linspace(.1,100000, 1000000)
In [384]:
log1 = np.log(x - x.min() +1) /np.log(x.max() - x.min() +1)
In [385]:
log2 = (np.log(x) - np.log(x.min()))/(np.log(x.max()) - np.log(x.min()))
In [392]:
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escalas logarítmicas normalizadas')
ax[0].plot(x, log2, 'b')
ax[0].set_title(r"$\frac{\log(x) - \log(x_{min})} {\log(x_{max}) - \log(x_{min})}$" )
ax[1].plot(x, log1, 'r')
ax[1].set_title(r"$\frac{\log(x - x_{min} +1)}{ \log(x_{max}- x_{min} +1)}$")
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
escala horizontal, ano sim, ano não¶
In [403]:
fig, ax = plt.subplots(2,1, sharey=True,sharex=True )
labels = years[::2]
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
plt.xticks(labels,labels)
fig.tight_layout()
plt.show()
Alternativa¶
In [411]:
fig, ax = plt.subplots(2,1, sharey=True)
labels = years[::2]
#plt.setp(ax, xticks=labels, xticklabels=labels)
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[0].set_xticks(labels)
ax[0].set_xticklabels(labels)
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
ax[1].set_xticks(labels)
ax[1].set_xticklabels(labels)
fig.tight_layout()
plt.show()
In [413]:
from scipy import stats
In [432]:
years_float = np.array([float(yr) for yr in years])
In [443]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels)
ax[i][j].set_title(biome_names_dict[i*3 + j])
ax[i][j].legend()
fig.tight_layout()
plt.show()
In [438]:
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html
Test seaborne style¶
In [439]:
import seaborn as sns; sns.set()
In [440]:
fig, ax = plt.subplots(2,1, sharey=True)
labels = years[::2]
#plt.setp(ax, xticks=labels, xticklabels=labels)
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[0].set_xticks(labels)
ax[0].set_xticklabels(labels)
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
ax[1].set_xticks(labels)
ax[1].set_xticklabels(labels)
fig.tight_layout()
plt.show()
In [441]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label="r-squared: {} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels)
ax[i][j].set_title(biome_names_dict[i*3 + j])
ax[i][j].legend()
fig.tight_layout()
plt.show()
In [653]:
%%HTML
<h1 id="1">PLOTS<h1>
Variação inter-anual. Totais.¶
In [509]:
fig, ax = plt.subplots()
labels = years[::2]
ax.plot(years, annual_total)
ax.set_xticks(labels)
ax.set_xticklabels(labels, fontsize=15)
ax.set_title('Interannual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
#plt.setp(labels_x,fontsize=15, family='Times New Roman')
plt.setp(labels_y,fontsize=15)
#plt.xlabel(r"Burned Pixels", fontsize = 12)
plt.ylabel(r"Burned Pixels", fontsize = 15)
plt.tight_layout()
plt.show()
In [640]:
#help(off)
Variação intra-anual. Totais¶
In [511]:
months = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split(' ')
In [513]:
fig, ax = plt.subplots()
ax.plot( np.arange(12),monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels", fontsize = 15)
plt.tight_layout()
plt.show()
In [520]:
residuals_monthly_total = (np.array(monthly_total) - mean(monthly_total))/mean(monthly_total)
In [539]:
fig, ax = plt.subplots()
ax.plot( np.arange(12),residuals_monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"$\frac{Burned Pixels - \langle Burned Pixels \rangle}{\langle Burned Pixels \rangle}$", fontsize = 15)
plt.tight_layout()
plt.show()
In [540]:
130070.096/3600
Out[540]:
$\frac{Burned Pixels - \langle Burned Pixels \rangle}{\langle Burned Pixels \rangle}$
Variação intermensal, totais¶
In [569]:
ym = itertools.product(years,months)
ymlabels = [ m + ' ' + y[-2:] for y, m in ym]
In [570]:
fig, ax = plt.subplots(figsize=(16,12))
labels = ymlabels[::6]
ax.plot( np.arange(12*17),all_months)
ax.axhline(y=mean(all_months), color='r', linestyle='--', label='average')
ax.set_xticks(np.arange(0,12*17,6))
ax.set_xticklabels(labels,fontsize=12, rotation=30)
labels_y = ax.get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax.set_title('Intra-monthly Variability',fontsize=15)
plt.legend()
plt.show()
Variação inter_anual em cada bioma¶
In [573]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
In [742]:
%%HTML
<h1 id="varintrabio"> Variação intra_anual em cada bioma</h1>
In [586]:
len(monthly_biome[0])
Out[586]:
In [587]:
ar_monthly_biome = np.array(monthly_biome)
In [589]:
ar_monthly_biome.shape
Out[589]:
In [595]:
monbio = ar_monthly_biome.T
In [598]:
monbio.shape
Out[598]:
In [600]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = months
for i in range(5):
for j in range(3):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(months, monbio[i*3 + j], 'o')
ax[i][j].set_xticks(np.arange(12))
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
In [601]:
fig, ax = plt.subplots(5,3, figsize =(18,12) )
labels = months
for i in range(5):
for j in range(3):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(months, monbio[i*3 + j], 'o')
ax[i][j].set_xticks(np.arange(12))
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
In [609]:
type(df_finall)
Out[609]:
In [657]:
biome_grouped.sum()/areas_biome
Out[657]:
In [639]:
biome_grouped.sum()[14][16*12]/areas_biome[14]
Out[639]:
In [658]:
biomes_pa = biome_grouped.sum()/areas_biome
In [681]:
biomes_pa.index.get_level_values(level=1)
Out[681]:
In [685]:
biomes_pa.index.get_level_values(level=1)
Out[685]:
In [690]:
intra_biome = biomes_pa.groupby(biomes_pa.index.get_level_values(level=1), sort=False).sum()
In [688]:
for a, b in biomes_pa.groupby(biomes_pa.index.get_level_values(level=1),sort=False):
print(a)
print(b)
In [691]:
intra_biome
Out[691]:
In [704]:
intra_biome[1]
Out[704]:
De novo¶
In [708]:
intra_biome = [ b for a, b in biomes_pa.groupby(biomes_pa.index.get_level_values(level=1),sort=False)]
In [709]:
intra_biome
Out[709]:
In [740]:
%%HTML
<h1 id="varintermonbio">Variação inter-anual por mês em cada bioma</h1>
In [738]:
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [695]:
calendar.month_name[3]
Out[695]:
In [745]:
sns.set_context("paper")
In [746]:
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [747]:
sns.set_style("whitegrid")
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [784]:
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
#slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
#ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
sns.regplot(years_float, data, ax=ax[i][j])
#ax[i][j].set_xticks(labels)
#ax[i][j].set_xticklabels(labels,fontsize=15)
#labels_y = ax[i][j].get_yticklabels()
#plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].set_ylabel('')
#ax[i][j].legend(fontsize=15)
#fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [765]:
x = np.linspace(1,200,10)
y = 4*x
df = pd.DataFrame(y, index=x)
In [766]:
df
Out[766]:
In [768]:
df.columns = ['y']
In [769]:
df.index.name= 'x'
In [771]:
df
Out[771]:
In [777]:
df =pd.DataFrame(dict(zip(['x','y'],[x,y])))
In [778]:
df
Out[778]:
In [783]:
ax= sns.regplot( x = 'x', y='y',data=df);ax.set(ylabel=''); plt.show()
In [796]:
pwd
Out[796]:
In [795]:
cd ../MODIS/MCD64A1
In [803]:
import datetime
from pyhdf.SD import SD, SDC
#import matplotlib.path as mplPath
import subprocess
import re
DATADIR = 'C:/Users/alpha/Documents/dados/MCD64A1/'
NECOS = 847
def get_tilehv(fn):
pattern = 'MCD64A1\.A\d{4}\d{3}\.h(\d{2})v(\d{2}).{22}'
# MCD64A1.A2005032.h14v04.006.2017017093958.hdf
matches= re.findall(pattern, fn)
if matches:
h, v = matches[0]
return int(h), int(v)
def get_coord(h,v):
lon_f = 'lon_coordh{}v{}.npy'.format(h,v)
lat_f = 'lat_coordh{}v{}.npy'.format(h,v)
lon = np.load(lon_f)
lat = np.load(lat_f)
return lon, lat
def get_julian(yr, mo):
'''
Return the first julian day of given month and year.
'''
dt = datetime.datetime(yr, mo, 1)
tt = dt.timetuple()
return tt.tm_yday
def julian2month(yr,jul):
# month from 1 to 12
values = [get_julian(yr,i) for i in range(1,13)]
#return month from 0 to 11
idx = [i for i in range(12) if values[i] == jul]
return idx[0]
def get_month(fn):
pattern = 'MCD64A1\.A(\d{4})(\d{3}).{29}'
matches= re.findall(pattern, fn)
year, julian = matches[0]
year, julian = int(year), int(julian)
return julian2month(year,julian)
file_ = DATADIR + 'MCD64A1.A2001152.h17v05.006.2017012103748.hdf'
h , v = get_tilehv(file_)
print('h: {}; v: {}'.format(h,v))
lon, lat = get_coord(h,v)
fich = SD(file_, SDC.READ)
# select sds
sds_obj = fich.select('Burn Date')
# get sds data
data = sds_obj.get()
attrs = sds_obj.attributes(full=1)
vra=attrs["valid_range"]
valid_range = vra[0]
fva=attrs["_FillValue"]
_FillValue = fva[0]
lat = lat.reshape(data.shape)
lon = lon.reshape(data.shape)
# Mask the data. logical_or from numpy
data = data.astype('double')
#invalid = logical_or(data < valid_range[0], data > valid_range[1])
invalid = np.logical_or(data < 1, data > valid_range[1])
invalid = np.logical_or(invalid, data == _FillValue)
# nan from numpy
data[invalid] = np.nan
# ma, isnan from numpy
data = np.ma.masked_array(data, np.isnan(data))
idx =np.logical_not(data.mask)
lat_val =lat[idx]
lon_val = lon[idx]
data_val = data[idx]
result = data_val.size
print('Result: {}'.format(result))
In [804]:
import pickle
from mpl_toolkits.basemap import Basemap
import matplotlib.path as mplPath
In [805]:
m = Basemap(projection='cyl')
shpf = m.readshapefile('ecoregions2017/ecoregions2017','ecos')
rings =[mplPath.Path(eco,closed=True) for eco in m.ecos]
ecos_info =m.ecos_info
In [809]:
result = np.zeros((NECOS),dtype=int)
In [807]:
br_vh = []
with open('bounded_rings.txt', 'rb') as f:
br_string = f.read()
br = pickle.loads(br_string)
br_vh = br[v,h]
In [810]:
not_included = []
for i, px in enumerate(data_val):
pt = ( lon_val[i], lat_val[i])
included = False
for ring_idx in br_vh:
ring = rings[ring_idx]
if ring.contains_point(pt):
ecoid = int(ecos_info[ring_idx]['ECO_ID'])
result[ecoid] += 1
included = True
break
if not included:
not_included.append(pt)
In [812]:
result.sum()
Out[812]:
In [813]:
len(not_included)
Out[813]:
Total por Ecoregiões¶
In [836]:
total_ecoreg = df_finall.iloc[:12*17].sum() / df_finall.iloc[12*17]
In [838]:
total_ecoreg
Out[838]:
In [839]:
ecos = m.ecos
In [859]:
cd -
In [860]:
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(total_ecoreg.values.astype('float'))/np.log1p(total_ecoreg.max())
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['ECO_ID']] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Index Burned Area 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('totalecologstd', dpi=1200, bbox_inches='tight')
In [861]:
#total_ecoreg.values.astype('float')
In [862]:
#np.log1p(total_ecoreg.values.astype('float'))
In [863]:
#total_ecoreg.values
In [858]:
pwd
Out[858]:
In [864]:
np.log1p(total_ecoreg.max())
Out[864]:
In [865]:
total_ecoreg.max()
Out[865]:
In [868]:
total_ecoreg.size
Out[868]:
In [883]:
areas_biome
Out[883]:
In [892]:
df_finall.iloc[:12*17].groupby(df_finall.loc['Atributos'].loc['BIOME_NUM'], axis=1).sum().sum()/areas_biome
Out[892]:
In [893]:
total_biome_area
Out[893]:
In [894]:
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = int(total_biome_area.max()) +1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
#total_trans = np.log(total_pa + 1)/np.log(total_pa.max() + 1)
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_biome_area[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_biome_area[0] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Burned Px / Area 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomesnotransform2', dpi=1200, bbox_inches='tight')
In [919]:
df_finall[:][:12*17].T
Out[919]:
In [921]:
df_finall.T['Atributos']['BIOME_NUM']
Out[921]:
In [928]:
lala = pd.concat([df_finall[:][:12*17].T,df_finall.T['Atributos']['BIOME_NUM']], axis=1)
In [951]:
start = datetime.datetime.now()
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(total_biome_area)/np.log1p(total_biome_area.max())
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform2', dpi=1200, bbox_inches='tight')
print('It took {} minutes.'.format((datetime.datetime.now()-start)/60))
In [950]:
total_trans
Out[950]:
In [952]:
import time
In [953]:
start = datetime.datetime.now()
time.sleep(29)
datetime.datetime.now()-start
Out[953]:
In [960]:
print(str(datetime.timedelta(0, 29, 10693)).format('%m:%s'))
In [968]:
datetime.timedelta(0, 29, 10693).seconds
Out[968]:
In [971]:
len([ df_finall[n*12: (n+1) *12].sum() for n in range(17)])
Out[971]:
In [973]:
n= 1
df_finall[n*12: (n+1) *12].sum()
Out[973]:
In [975]:
len(areas)
Out[975]:
In [977]:
intereco =[ df_finall[n*12: (n+1) *12].sum()/areas for n in range(17)]
In [992]:
intereco_max= max(map(max,intereco))
In [998]:
for yr in range(17):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(intereco[yr])/np.log1p(intereco_max)
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['ECO_ID']] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Index Burned Area ' + years[yr]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('ecologstd' + years[yr], dpi=1200, bbox_inches='tight')
In [1028]:
annual_biome
Out[1028]:
In [1035]:
type(annual_biome[0])
Out[1035]:
In [1054]:
denominador = np.log1p(max(map(max,annual_biome)))
for yr in range(17):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(annual_biome[yr])/denominador
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.loc[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.iloc[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned' + years[yr]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform' + years[yr], dpi=1200, bbox_inches='tight')
In [1053]:
annual_biome[0]
Out[1053]:
In [1039]:
max(map(max,annual_biome))
Out[1039]:
In [1043]:
open('df_finall','wb').write(pickle.dumps(df_finall))
Out[1043]:
In [1044]:
lala = pickle.loads(open('df_finall','rb').read())
In [1047]:
pwd
Out[1047]:
In [1051]:
int(ecos_info[0]['BIOME_NUM'])
Out[1051]:
In [1055]:
monthly_biome
Out[1055]:
In [1061]:
notincluded2013 = pickle.loads(open('notincluded2013', 'rb').read())
In [1062]:
notincluded2014 = pickle.loads(open('notincluded2014', 'rb').read())
In [1065]:
len(notincluded2014)
Out[1065]:
In [1068]:
lons, lats = zip(*notincluded2013)
In [1078]:
lons=list(lons)
lats= list (lats)
In [1084]:
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
x, y = map(lons, lats)
map.plot(lons,lats, 'mo', markersize=5)
plt.show()
In [1077]:
lons, lats = zip(*notincluded2014)
Out[1077]:
In [1085]:
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
x, y = map(lons, lats)
map.plot(lons,lats, 'mo', markersize=5)
plt.show()
In [1095]:
lons14, lats14 = zip(*notincluded2014)
lons13, lats13 = zip(*notincluded2013)
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.readshapefile('ecoregions2017/ecoregions2017','ecos',linewidth=0.1)
x, y = map(lons13, lats13)
map.plot(lons13,lats13, 'ro', markersize=2, label='2013')
map.plot(lons14,lats14, 'b+', markersize=2, label='2014')
plt.legend()
plt.show()
In [1103]:
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='white',lake_color='aqua')
map.drawcoastlines()
map.plot(lons13,lats13, 'go', markersize=2, label='2013')
map.plot(lons14,lats14, 'ro', markersize=2, label='2014')
plt.legend()
plt.show()
In [1091]:
monthly_biome
Out[1091]:
In [1099]:
max([ mb.max() for mb in monthly_biome])
Out[1099]:
In [1100]:
np.log1p(20120.9686490036)
Out[1100]:
In [1092]:
denominador2 = 9.9095674648050416
In [1093]:
for mo in range(12):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(monthly_biome[mo])/denominador2
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.loc[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.iloc[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned ' + calendar.month_name[mo + 1]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform' + calendar.month_name[mo + 1], dpi=1200, bbox_inches='tight')
